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DEVELOPMENT  OF  AN  ADAPTIVE  BOUNDARY-FITTED  COORDINATE  CODE 
FOR  USE  IN  COASTAL  AND  ESTUARINE  AREAS 

i 

PART  I:  INTRODUCTION  ■ 

i 

I 

i 

1.  The  mathematical  modeling  of  the  hydrodynamics  of  a  body  of  water 
plus  the  transport  and  dispersion  of  a  constituent  within  that  body  involves 

the  solution  of  a  set  of  partial  differential  equations  expressing  the  conser-  ] 

vation  of  mass,  momentum,  and  energy  of  the  flow  field  along  with  a  transport 
equation  for  the  constituent.  These  equations  involve  derivatives  with  re¬ 
spect  to  time  as  well  as  three  spatial  dimensions.  However,  a  simplification 
that  is  often  made  in  treating  relatively  shallow  bodies  of  water  which  are 
well  mixed  over  the  depth  is  to  vertically  average  the  three-dimensional  (3~D) 
equations  to  yield  a  two-dimensional  (2-D)  set  for  nearly  horizontal  flows. 

Numerical  Techniques 

2.  Since  the  governing  equations  are  nonlinear,  analytic  solutions  in 
general  cannot  be  found  and  one  is  forced  to  resort  to  numerical  techniques  to 
obtain  solutions.  The  two  most  common  such  techniques  are  the  finite  differ¬ 
ence  method  ( FDM)  and  the  finite  element  method  (FEM).  There  are,  of  course, 
both  advantages  and  disadvantages  to  each  of  these  approaches. 

3.  Perhaps  the  most  often  quoted  advantage  of  the  FEM  is  that  with  this 
approach,  physical  boundaries  coincide  with  computational  net  points.  There¬ 
fore  the  modeling  of  flow  within  an  irregular  domain  can  be  more  accurately 
handled  than  with  the  older  rectangular  FDM  wherein  the  approach  is  to  con¬ 
struct  a  rectangular  grid  over  the  domain,  which  forces  the  boundaries  to  be 
represented  in  a  "stair-stepped"  fashion.  However,  a  disadvantage  of  many 
finite  element  models  is  the  excessive  computational  time  required  compared 
with  typical  finite  difference  models  having  the  same  number  of  mesh  points. 

As  discussed  by  Johnson,  Thompson,  and  Baker  (198^)  this  occurs  because  of  the 
manner  in  which  the  system  of  resulting  algebraic  equations  is  usually  solve! 
in  finite  element  models.  An  additional  disadvantange  is  that  the  FEM  is  more 
cumbersome  to  code  into  a  computer  model  than  the  FDM.  This  can  be  a  problem 
not  only  during  the  development  of  the  model  but  also  can  inereise  the  level 
of  effort  required  during  later  model  modifications. 
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Boundary-Fitted  Coordinates  Concept 


4.  Accepting  that  the  FDM  possesses  an  advantage  in  simplicity  and  per¬ 
haps  computational  costs,  a  logical  question  is  whether  or  not  one  can  develop 
ways  to  circumvent  the  major  disadvantage  of  having  to  represent  irregular 
boundaries  in  a  staii — stepped  fashion.  One  method  is  to  make  computations  on 
curvilinear  grids  so  that  one  of  the  curvilinear  coordinates  always  follows  a 
boundary.  Vertically  averaged  hydrodynamic  models  developed  by  Johnson  (1980) 
and  Sheng  and  Hirsh  (1984)  are  examples  of  models  employing  a  general  non- 
orthogonal  curvilinear  grid.  Through  coordinate  transf ormations ,  irregular 
boundaries  and  variable  grid  spacing  can  be  more  accurately  handled  while 
still  making  use  of  the  simplicity  of  finite  differences  to  obtain  solutions. 
Since  the  boundary-fitted  coordinate  system  has  a  coordinate  line  coincident 
with  all  boundaries,  all  boundary  conditions  can  be  expressed  at  grid  points; 
and  normal  derivatives  can  be  represented  using  only  finite  differences  be¬ 
tween  grid  points  on  coordinate  lines.  No  interpolation  is  needed,  even 
though  the  coordinate  system  is  not  orthogonal  at  the  boundary.  A  general 
discussion  is  given  by  Thompson,  Warsi  ,  and  Mastin  (1985). 

5.  To  enable  efficient  flow-transport  modeling  on  boundary-fitted 
grids,  a  numerical  generator  is  needed  to  provide  the  (x,y)  location  of  the 
curvilinear  coordinates.  One  of  the  earliest  grid  generators  W3s  developed  by 
Thompson,  Thames,  and  Mastin  (1977)  with  a  later  version  called  WESCOR  pro¬ 
vided  by  Thompson  H  983). 

WESCOR 

b.  The  WESCOR  code  generates  a  boundary-conforming  curvilinear  coor- 
linste  system  for  a  general  J-D  region  with  boundaries  of  arbitrary  shape  and 
with  boundary  intrusions  and  internal  obstacles,  such  as  islands,  arbitrary  in 
shape  and  number.  The  grid  is  generated  from  the  numerical  solution  of  a  set 
of  elliptic  partial  differential  equations  by  accelerated  point  Gauss-Seidel 
iteration. 

7 .  These  equations  are  written  in  the  transformed  space,  which  is  in¬ 
herently  "*'ot angular  with  a  •.quar'e  grid.  All  computations,  both  to  generate 
tiie  g  ■  i  1  and  subsequently  to  so!  v*>  partial  differential  equation:;  for  physical 
problems  on  the  grid,  ar*»  lone  in  this  transformed  space,  so  that  all  boundary 


•v,v. 


extension  of  the  control  functions  in  the  original  elliptic  generation  in 
WESCOR. 

11.  With  either  adaptive  mechanism  the  code  first  generates  a  grid 
using  the  original  WESCOR  system  without  regard  to  depths.  The  input  depth 
distribution,  which  is  specified  by  an  arbitrary  arrangement  of  points  with  no 
order  or  pattern,  is  then  interpolated  onto  this  initial  grid.  No  further 
interpolation  is  needed  as  the  adaptive  grid  generation  system  is  solved  for  a 
grid  that  is  concentrated  according  to  increasing  depth. 

12.  In  the  following  sections,  discussion  of  the  theory  of  the  adaption 
mechanism,  the  interpolation  procedures,  input  required,  and  finally  an  exam¬ 
ple  application  are  presented.  User  instructions  and  sample  job  streams  are 
presented  in  the  appendices. 


PART  II:  ADAPTIVE  GRIDS 


13.  The  generation  of  adaptive  grids  is  discussed  in  detail  by 
Thompson,  Warsi,  and  Mastin  (1985)  and  a  recent  survey  is  given  by  Thompson 
(1985).  A  discussion  of  potential  applications  in  numerical  hydrodynamic 
modeling  is  given  by  Johnson,  Thompson,  and  Baker  (1984).  Briefly,  an  adap¬ 
tive  grid  senses  gradients  of  some  physical  quantity  and  adjusts  itself  auto¬ 
matically  to  better  resolve  those  gradients,  i.e.,  to  concentrate  lines  in  the 
high  gradient  regions.  (The  adaption  here  is  through  movement  of  the  grid 
points,  not  the  addition  of  more  points.) 

14.  With  the  time  derivatives  at  fixed  values  of  the  physical  coordi¬ 
nates  transformed  to  time  derivatives  taken  at  fixed  values  of  the  curvilinear 
coordinates,  no  interpolation  is  required  when  the  adaptive  grid  moves.  The 
time  derivative  transforms  as  follows: 


(1  ) 


where  x,y  are  the  components  of  the  grid  speed.  The  computation  thus  can  be 
done  on  a  fixed  grid  in  the  transformed  space,  without  need  of  interpolation, 
even  though  the  grid  points  are  in  motion  in  physical  space.  The  influence  of 
the  motion  of  the  grid  points  is  registered  through  the  grid  speed  appearing 
in  the  transformed  time  derivative. 

15.  Several  considerations  are  involved  here,  some  of  which  are  con¬ 
flicting.  The  points  must  concentrate,  and  yet  no  region  can  be  allowed  to 
become  devoid  of  points.  The  distribution  also  must  retain  a  sufficient 
degree  of  smoothness,  and  the  grid  must  not  become  too  skewed,  else  the  trun¬ 
cation  error’  of  equations  solved  on  the  grid  will  be  increased.  This  means 
that  points  must  not  move  independently,  but  rather  each  point  must  somehow  be 
coupled  at  least  to  its  neighbors.  Also,  the  grid  points  must  not  move  too 
far  or  too  fast,  else  oscillations  may  occur. 


Variational  Formulation 


16.  Thus,  on  the  one  hand,  there  is  a  need  to  force  grid  points  to  con¬ 


centrate  near  large  solution  gradients;  on  the  other  hand,  there  is  a  need  to 
generate  grids  that  are  relatively  smooth  and  do  not  deviate  too  much  from 


being  orthogonal.  The  calculus  of  variations  is  well  suited  to  handle  such 
problems  since  integrals  over  the  grid  can  be  written  that  measure  the  three 
desired  features  discussed  above,  namely,  (a)  a  concentration  of  grid  points 
near  large  solution  gradients,  (b)  a  smooth  distribution  of  grid  points,  and 
(c)  a  relatively  orthogonal  grid.  The  variational  formulation  of  adaptive 
grids  is  discussed  in  detail  by  Brackbill  and  Saltzman  (1982)  and  by  Thompson, 
Warsi,  and  Mastin  (1985).  A  brief  discussion  follows. 

17.  The  area  of  a  computational  cell  in  two  dimensions  (the  volume  in 
three  dimensions)  is  given  by  the  Jacobian,  J  ,  of  the  mapping: 


J  ■  Vn  ■  Vc 


Therefore,  if  the  integral 


I  =  fw(x,y)  J  dA  ,  (3) 

W  D 

which  is  a  measure  of  the  weighted  variation  of  the  computational  cell  size 
over  the  grid,  is  minimized  for  some  weight  function  w(x,y)  ,  which  is 
related  to  the  solution  gradient,  a  concentration  of  grid  points  in  high 
gradient  regions  can  be  obtained.  In  other  words,  where  the  weight  function 
w(x,y)  becomes  large  the  cell  size  becomes  small,  and  where  w(x,y)  becomes 
small  the  size  of  the  computational  cell  becomes  large. 

18.  Likewise,  the  smoothness  of  the  grid  is  measured  by  the  integral 

I  =  f  [m)2  +  ( Vn)  |dA  (4) 

3  DL  J 

with  the  orthogonality  measured  by 


I 

o 


J(V£  •  Vn)2  J3  dA 
D 


(5) 


where  the 


j3 


is  added  to  cause  orthogonal i ty  to  be  emphasized  more  in  larger 


cells.  Note  that  V£  •  Vn  =  0  for  an  orthogonal  grid.  Therefore  an  adaptive 
grid  generator  can  be  developed  by  minimizing  a  sum  of  the  three  integrals 
given  in  Equations  3~5,  i.e., 


19.  From  the  above  expressions  for  T  ,  T  ,  and  I  we  have  the 

3  0  W 

dimensional  relations 
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where 

N  =  character istic  number  of  points 

L  =  character istic  length 

W  =  average  weight  function  over  the  field 


W  =  j  jw(x,y)  dA 
A  D 

Thus,  for  Equation  6  to  be  dimensionally  correct,  we  take 


(7) 


where  A'  and  A'  are  positive  constants.  The  character ist io  length  and 
o  w 

number  of  points  might  logically  be  taken  as  the  square  roots  of  the  total 
area  and  the  total  number  of  points  in  the  field,  respectively. 

20.  Obviously,  by  selecting  appropriate  values  of  A^  ,  A^  ,  and  A^  , 
grids  that  emphasize  smoothness,  concentration  of  grid  points,  or  orthog¬ 
onality  can  be  generated.  The  variational  formulation  thus  provides  a  com¬ 
petitive  enhancement  of  these  three  grid  character istics .  A  note  of  caution 
concerning  orthogonality  is  perhaps  needed.  Purely  orthogonal  grids  cannot  be 
generated  when  prescribing  the  location  of  boundary  points,  even  if  the  con¬ 
dition  •  Vn  =  0  is  enforced.  Since  derivatives  of  the  coordinates  have 
to  be  specified  to  satisfy  orthogonal i ty  at  the  boundaries,  specifying  the 
location  of  the  points  overspecifies  the  problem  with  second-order  partial 
differential  equations  as  the  grid  generation  system.  Therefore,  in  order  to 
generate  a  strictly  orthogonal  grid,  the  boundary  points  must  be  t L lowed  to 
move  on  the  boundary.  In  many  cases  this  constitutes  a  major  restriction  on 

the  generation  of  a  useful  grid.  Therefore  large  values  of  compared  wit.!: 

o 

values  of  A'  and  A'  will  only  enhance  tne  orthogonal i ty  of  tne  grid.  In 

3  W 

general  the  grid  is  still  nonorthogonal ,  and  all  terms  in  the  governing  trans¬ 


formed  partial  differential  equations  must  be  retained. 

21.  For  a  2-D  adaptive  grid  generator,  partial  differential  equitions 


for  the  solution  for  the  physical  coordinates  (x,y)  are  desired.  These  are 
obtained  through  the  calculus  of  variations  applied  to  the  minimization  of  the 
sum  of  integrals  in  Equation  6.  In  general ,  if  we  wish  to  find  functions 
d> ^ ( x , y )  that  are  differentiable  on  (x,y)  and  satisfy  fixed  constraints  on 
the  boundary  of  the  domain  that  minimize  some  integral  functional 


f  (  3*i  5M 

Jf  W’V  st'  3 r/dx  dy 


the  calculus  of  variations  gives  that  <?>.(x,y)  are  the  solutions  of  the 
Euler-Lagrange  equations 


r  af  1  _  3f_ 

^ 3 ( V<}> ^  ) J  9^ 


22.  The  2-D  Euler-Lagrange  equations  for  the  minimization  of  the  sum  of 
integrals  in  Equation  6  are  given  by  Brackbill  and  Saltzman  (1982)  as 
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For  completeness,  the  coefficients  are  reproduced  from  the  cited  paper  and 
presented  below: 
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Equation  10  constitutes  the  variational  adaptive  grid  generator  used  in  the 
present  code. 

Control  Function  Formulation 

23.  The  elliptic  generation  system  used  in  the  WESCOR  code  is  as 
follows: 
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where  now  the  coefficients  are  redefined  as 
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Here  P  and  Q  are  the  control  functions  that  serve  to  concentrate  the  grid 
lines . 


2H.  In  one  dimension,  with  =  0  ,  these  equations  reduce  to 
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Now  if  the  product  of  a  weight  function  w  and  the  grid  spacing  is  equally 
distributed  over  the  grid  we  have,  in  one  dimension, 
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Then  by  analogy  we  make  the  connection 


Generalizing  we  have 


(12 


Equation  11,  with  the  control  functions  evaluated  from  Equation  12,  consti¬ 
tutes  the  adaptive  generation  system  based  on  the  original  elliptic  system  of 
the  WESCOR  code. 


PART  III:  DEPTH  INTERPOLATION 


25.  A  basic  input  to  depth-averaged  hydrodynamic  models  is  the  water 
depth  associated  with  each  computational  cell.  These  depths,  relative  to  some 
datum,  are  normally  obtained  from  National  Oceanic  and  Atmospheric  Administra¬ 
tion  ( NOAA )  charts.  An  efficient  grid  generator  that  forces  grid  points  to 
cluster  in  navigation  channels  through  adaption  to  the  water  depth  requires 
the  interpolation  of  the  input  depth  field  determined  from  the  NOAA  charts. 

The  number  of  input  depth  points  and  their  distribution  should  be  arbitrary. 
Various  interpolation  schemes  that  have  been  incorporated  are  discussed  below. 


Taylor  Series  Expansion  About  Nearest  Depth  Point 


26.  With  (x,y)  the  nearest  depth  point  to  a  grid  point  (x,y),  the  depth 
f(x,y)  at  the  grid  point  is  given  by  a  Taylor  series  expansion  about  the 
depth  point  as  follows: 

f(x,y)  =  f(x,y)  +  fx(x  -  x)  +  fy(y  -  y)  «•  j  fxx(x  -  x)2 

+  \  fyy(y  '  y)2  +  fxy(x  '  x)  {y  ~  y)  +  °3) 

where  all  derivatives  are  evaluated  at  the  depth  point. 

I  27.  These  derivatives  are  determined  by  Taylor  series  expansions  of  the 

given  depths  at  depth  points  neighboring  the  one  in  question  about  the  latter 

(cf.  Figure  1).  Thus  if  the  five  nearest  depth  points  to  the  depth  point 
(x,y)  are  (x^yj)  i  =  1,2,3,^»5  ,  we  then  have 
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Since  the  left  sides  here  are  known  as  values  at  depth  points,  this  system 
constitutes  five  equations  for  the  five  derivatives  at  the  depth  point  (x,y). 
If  fewer  derivatives  are  included  in  the  expansion  about  the  depth  point  then 


the  number  of  equations  in  the  system  (I1*)  is  reduced  accordingly.  For 

example,  for  a  linear  expansion,  only  f  and  f  are  used  in  Equation  13, 

*  y 

so  only  two  neighboring  depth  points  are  used  and  the  system  (14)  consists  of 

two  equations  only.  The  code  provides  for  the  use  of  two  (f  ,f  ),  three 

x  y 

(fx,fy,fxy),  four  (fx,f  ,fxx,fyy) ,  and  all  five  derivatives. 

28.  This  interpolation  is  implemented  by  first  sweeping  all  the  depth 
points,  at  each  of  which  the  appropriate  numbering  of  neighboring  points  is 
located  and  the  required  derivatives  are  evaluated  by  the  solution  of  the 
system  (14).  Then  the  grid  points  are  swept,  determining  the  nearest  depth 
point  at  each  and  determining  the  depth  from  Equation  13. 

Interpolation  Within  Triangles  or  Quadr i laterals  of  Depth  Points 

29.  If  the  depth  at  a  grid  point  (x,y)  is  written  as  the  linear 
f  unction 


f(x,y)  =  a  +  bx  +  cy  (15) 

then  the  coefficients  a  ,  b  ,  and  c  can  be  determined  by  evaluating  Equa¬ 
tion  15  at  three  depth  points  surrounding  the  grid  point  in  question  (cf. 
Figure  2).  Thus  we  have  the  system 

f  i  =  a  +  bXj:  +  cyt  i  =  1 ,2,3  (16) 

of  three  equations  for  the  three  coefficients. 

30.  Similarly,  with  the  depth  given  by  the  product  of  linear  functions 
in  each  coordinate  we  have 

f ( x , y )  =  (a  *  bx)  (c  +  dy) 

=  ac  +  bex  +  ady  +  bdxy 

Redefining  the  coefficients,  this  can  be  written  as 

f(x,y)  =  a  *  bx  *  cy  +  dxy  (17) 

and  now  evaluation  at  four  surrounding  points  gives  the  system 


=  a  +  +  cy^^  +  dx^ 
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(18) 


31.  Implementation  proceeds  by  sweeping  the  grid  points,  at  each  of 
wMch  the  three  or  four  nearest  surrounding  depth  points  are  located  and  the 
coefficients  are  determined  by  solution  of  the  system  (16)  or  (18).  Then  the 
depth  at  the  grid  point  is  evaluated  from  Equation  15  or  Equation  17. 


Inverse-Power  Interpolation  Among  Depth  Points 


32,  Here  the  depth  at  a  grid  point  (x,y)  is  written  in  terms  of  inverse 
powers  of  the  distance  to  each  point  of  a  surrounding  group  of  nearest  depth 
points  (cf.  Figure  3)  as  follows: 
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where  d^  =  -  xj2  +  (y  -  y.j2  *3  the  distance  from  the  grid  point  in 

question  to  the  depth  point  (x^.y^)  and  m  is  a  specified  power,  one  or 
greater.  The  code  provides  for  two,  three,  four,  five,  or  all  depth  points  to 
be  used  in  Equation  19.  The  exponent  m  is  an  input  quantity,  higher  values 
of  this  exponent  giving  sharper  variations  in  depth. 


Location  of  Neighboring  Depth  Points 


33.  The  accuracy  of  the  interpolation  is  enhanced  if  the  group  of  depth 
points  used  generally  surrounds  the  point  in  question.  This  requires  a  some¬ 
what  different  selection  procedure  for  each  succeeding  point  chosen  as  de¬ 
scribed  below.  Each  selected  point  is,  of  course,  excluded  from  the  selection 
procedure  for  the  succeeding  points. 

First  point 


34.  The  first  depth  point  selected  is  simply  the  nearest  one  to  the 
point  in  question:  (The  notation  in  all  the  following  figures  is  that  the 
vector  r  points  from  point  i  to  point  j  ,  with  the  point  about  which 

the  group  is  being  formed  designated  as  point  0  .) 


unit  vector  a  along  the  bisector  of  the  angle  between  the  first  two  points 
is  given  by 
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Then  the  third  point  will  lie  in  the  shaded  region  if 
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which,  upon  substitution  for  a  ,  becomes 
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37.  The  nearest  point  satisfying  this  criterion  is,  of  course,  chosen. 
Selection  of  the  third  point  in  this  manner  assures  that  ‘he  point  about  which 
the  group  is  formed  will  lie  inside  a  triangle  formed  by  the  three  depth 
points  chosen. 

Fourth  point 

38.  The  fourth  point  selected  is  required  to  lie  outside  the  /.ones 
formed  by  extending  the  sides  of  the  triangle  of  points  already  chosen: 


In  the  diagram  here,  a  is  the  unit  vector  on  the  bisector  of  the  angle  be¬ 
tween  the  vector  from  point  1  to  points  2  and  3-  Thus 
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The  fourth  point  will  lie  in  the  region  between  the  extensions  of  the  vectors 
a12  and  ^  ,  as  shown,  if 
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Following  a  development  analogous  to  that  given  ahov“  in  connection  with  the 
third  point,  this  becomes 
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This,  however,  does  not  ensure  that  the  fourth  point  is  outside  the  triangle 
formed  by  the  first  three  points.  Therefore  we  require  in  addition  that 
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Crossing  a ^  into  both  sides  to  eliminate  q  ,  we  have 
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Point  4  then  will  lie  outside  the  triangle  if 
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The  inequalities  (22)  and  (23)  apply  for  a  point  outside  the  triangle  and  be¬ 
tween  the  extensions  of  the  sides  with  a  vertex  at  point  1.  Generalizing  to 
the  other  two  vertices,  the  fourth  point  is  selected  as  the  nearest  point 
satisfying  conditions  (24)  and  (25)  below  for  any  of  the  three  values  of  i  , 
i.e.,  for  any  vertex  of  the  triangle: 
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(i,m,n)  cyclic 


Selection  of  the  fourth  point  in  this  manner  guarantees  that  the  four  points 
form  a  convex  quadrilateral  about  the  point  in  question. 

Fifth  point 

39.  Here  again  the  nearest  point  is  chosen,  without  regard  to  other  con 
siderations.  This  simple  choice  of  the  fifth  point  does  not  necessarily  pro¬ 
duce  a  convex  polygon,  but  the  extra  complications  involved  in  doing  so  was  no 
considered  justified  in  view  of  the  lack  of  real  incentive  to  use  five  points. 


Additional  groups 

ilO.  In  some  cases  with  fewer  depth  points  it  is  possible  that  the 
interpolation  is  inaccurate  because  the  particular  group  of  points  selected 
gives  a  false  impression  of  very  large  gradients.  Therefore  provision  is  made 
for  selecting  more  than  one  group  of  points  and  basing  the  interpolation  on 
the  group  giving  the  smallest  gradients.  The  code  also  automatically  replaces 
points  in  any  group  leading  to  a  zero  determinant  in  the  solution  of  the 
linear  system  for  the  derivativ  s  or  coefficients. 


PART  IV:  GENERATION  PROCEDURE 


Iterative  Solution 


HI.  The  numerical  solution  of  Equation  10  or  Equation  11  is  accom¬ 
plished  in  finite  difference  form  by  replacing  all  derivatives  with  second- 
order  central  difference  expressions.  The  resulting  difference  equations  are 
nonlinear  since  the  coefficients  of  the  second  derivatives  involve  the  first 
derivatives.  This  system  of  difference  equations  is  solved  by  accelerated 
point  Gauss-Seidel  iteration.  In  both  the  first  and  second  derivatives,  all 
off-center  values  above,  or  to  the  right  of,  the  point  of  evaluation  are 
evaluated  at  the  previous  iteration,  while  those  below,  or  to  the  left,  are 
evaluated  at  the  present  iteration.  With  the  central  values  in  the  second 
derivatives  factored  together,  we  then  have  a  2*2  system  for  x  and  y  at 
the  point  in  the  case  of  Equation  10  and  an  uncoupled  system  with  Equation  11. 
The  field  is  swept  repetitively  until  convergence  to  some  prescribed  tolerance 

Weight  Function 


42.  The  components  of  the  gradient  of  the  weight  function  are  given  by 
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and  the  values  of  the  weight  function  at  each  point  are  continually  updated 
during  the  iteration  according  to  Equation  1.  Since  the  weight  function  is 
not  time-dependent  on  the  physical  field  in  the  present  application  we  have 


=  0 


so  that  Equation  1  reduces  to 


Then  using  first-order  forward  difference  expressions  for  the  time  derivatives 
and  canceling  the  At  we  have 


Aw  =  w  Ax  ♦  w  Ay 
x  y 


(27 


which  could  also  have  been  obtained,  of  course,  by  direct  application  of  the 

chain  rule.  Equation  27  then  serves  to  determine  the  change  in  w  at  each 

point  in  terms  of  the  change  in  x  and  y  at  each  iteration. 

^3.  The  code  first  generates  a  grid  from  the  original  WESCOR  system 
without  adaption.  Depth  values  are  then  interpolated  onto  this  initial  grid 
from  the  set  of  depth  points.  The  depth  points  may  be  input  in  any  order  and 
without  any  pattern.  The  code  automatically  adds  the  points  on  the  outer 

boundary  to  the  input  set  of  depth  points,  with  a  depth  of  zero  for  the  for¬ 

mer,  and  eliminates  duplicated  points.  It  is  this  augmented  set  that  is  used 
in  the  interpolation.  These  zero  depths  on  the  outer  boundary  can  be  over¬ 
ridden  by  including  the  outer  points  that  are  to  have  nonzero  depth  in  the 
input  set  of  depth  points.  Several  forms  of  interpolation  are  provided  as 
discussed  in  PART  [[I.  After  the  interpolation  the  depths  can  be  smoothed  on 
the  initial  grid  points. 

44.  The  smoothed  interpolated  values  of  the  depth  on  the  initial  grid 
then  are  used  as  the  weight  function  in  the  generation  system,  the  values  of 
the  depth  being  continually  updated  at  each  point  from  Equation  13  as  the  grid 
moves  without  further  interpolation. 


Input 


45.  The  depth-adaptive  code,  called  WESCORA,  has  the  same  structure  as 
the  earlier  WESCOR  code  capable  of  generating  2-D  grids  on  arbitrary  regions 
with  interior  obstacles.  The  input  is  the  same  as  that  described  for  the 
WESCOR  code  by  Thompson  (1983),  with  a  few  additions  described  below.  Com¬ 
plete  input  instructions  for  WESCORA  are  given  in  Appendix  A. 

46.  The  total  number  of  depth  points  to  be  read  in  for  grid  adaption  is 
input  as  NDEP,  and  a  flag  NDOP  controls  the  printing  of  the  depth  points.  (If 
NDEP=0  the  code  runs  the  nonadaptive  system  only.)  The  flag  NDEPF  causes  the 
depth  points  to  be  read  from  file  12  instead  of  from  the  input.  The  type  of 
Interpolation  is  specified  by  NINP  as  follows: 
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If  adjusted  depths  are  used  for  the  adaption,  then  an  additional  depth  file 
containing  actual  depths  is  input  to  provide  real  depths  at  the  center  of  each 
computational  cell.  The  flag  LDEPF  causes  the  second  set  of  depth  points  to 
be  read  from  file  16  instead  of  from  the  input,  with  LDEP  points  contained  in 
the  file.  The  interpolation  of  the  actual  depth  field  onto  the  final  adaptive 
grid  can  be  a  different  type  from  that  employed  in  the  grid  generation.  The 
types  available  are  the  same  as  those  described  above.  The  use  of  an  adjusted 
depth  field  is  discussed  in  more  detail  in  PART  V. 

47.  The  parameter  NCOM  specifies  the  number  of  groups  of  points  to  be 
considered  in  the  interpolation  (the  default  is  one  group).  The  parameter 
NTES  suppresses  the  requirement  that  the  points  selected  form  a  convex  polygon 
about  the  point  in  question,  and  the  parameter  NSMO  suppresses  the  smoothing 
of  the  depths  after  interpolation  onto  the  grid  points. 

48.  The  three  parameters  controlling  the  relative  emphasis  on  concen¬ 
tration,  smoothness,  and  orthogonality  are  input  as  WFACI,  SFAC,  and  OFACI, 
respectively,  the  nominal  ranges  being  0-1.  (A  negative  value  of  WFACI  will 
cause  the  control  function  adaption,  rather  than  the  variational  adaption,  to 
be  used,  with  the  magnitude  in  the  nominal  range.  In  this  case  the  other  two 
parameters  are  irrelevant.) 

49.  The  depth  points  used  in  the  adaption  are  read  in  either  from  the 
input  or  from  file  12,  immediately  after  the  boundary  points  are  read,  with 
the  Cartesian  coordinates  of  the  point  in  columns  0-10  and  11-20,  and  the 
depths  in  columns  21-30.  The  depth  points  may  be  in  any  order,  and  duplica¬ 
tions  will  be  automatically  eliminated  by  the  code.  If  a  second  set  of  depth 


points  is  input,  they  are  read  in  either  from  the  input  or  from  file  16.  The 
format  is  the  same  as  noted  above. 


50.  An  additional  feature  has  also  been  added  allowing  the  use  of 
Neumann  boundary  conditions  to  produce  orthogonality  on  straight  boundaries  if 
desired.  This  feature  is  activated  by  NEUBOD  for  slab  sides  and  by  NEUOUT  for 
the  outer  boundary. 
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PART  V:  CHESAPEAKE  BAY  APPLICATION 


51.  For  demonstration  purposes,  the  adaptive  grid  generation  code 
WESCORA  has  been  applied  to  Chesapeake  Bay  (Figure  4). 

Chesapeake  Bay 

52.  Chesapeake  Bay,  located  on  the  east  coast  of  the  United  States,  is 
one  of  the  largest  estuaries  in  the  world.  The  main  bay  extends  approximately 
190  miles  north  from  the  ocean  entrance  in  the  Commonwealth  of  Virginia,  be¬ 
tween  Cape  Henry  and  Cape  Charles,  to  the  Susquehanna  River  in  the  State  of 
Maryland.  The  average  depth  of  the  bay  is  about  28  ft,  although  a  natural 
channel  with  depths  greater  than  50  ft  traverses  the  bay  for  more  than  60  pei — 
cent  of  its  length.  The  maximum  depth  of  175  ft  is  located  in  the  upper  bay 
near  Bloody  Point,  Kent  Island,  Maryland. 

53.  Like  many  coastal  plain  estuaries,  the  bay  is  irregular  in  shape 
varying  in  width  from  4  miles,  between  Annapolis  and  Kent  Island,  to  30  miles, 
in  the  middle  bay  off  the  Potomac  River.  More  than  64,000  square  miles  of 
drainage  area  empty  into  the  bay  through  more  than  50  different  tributary 
systems.  Five  major  western  shore  rivers  (Susquehanna,  Potomac,  James,  York, 
and  Rappahannock)  provide  approximately  90  percent  of  the  annual  freshwater 
discharge . 


Selection  of  Boundary  Points 


54.  A  basic  input  to  the  grid  generation  code  is  the  specification  of 
the  (x,y)  coordinates  of  the  boundary  points  (Figure  5).  The  degree  of  reso¬ 
lution  of  boundary  features  will,  of  course,  depend  upon  the  number  of  6 
and  n  lines  selected  in  the  transformed  rectangular  plane  as  well  as  the 
location  of  the  boundary  points.  The  actual  boundary  points  selected  for  use 
in  WESCORA  are  listed  in  Appendix  B.  Figure  5  is  presented  only  for  demon¬ 
stration  purposes.  The  rectangular  (£,n)  plane  corresponding  to  the  actual 
boundary  points  listed  in  Appendix  B  is  shown  in  Figure  6. 


55.  Initially,  depth  points  were  randomly  selected,  with  approximately 
800  points  specified  to  cover  the  navigation  channel  as  well  as  the  nonchannel 
areas.  Depths  ranging  from  0  on  the  boundary  to  a  maximum  of  156  ft  were  in¬ 
put.  All  of  the  interpolation  schemes  previously  discussed  were  tried  with 
little  success  in  forcing  coordinate  lines  to  closely  follow  the  navigation 
channel.  The  grid  presented  in  Figure  7,  which  resulted  from  a  3-point  in¬ 
verse  power  interpolation  with  an  exponent  of  4  (Equation  19),  does  show  some 
adaption  to  the  channel.  However,  the  depth  field  interpolated  on  the  initial 
grid  presented  in  Figure  8  evidently  did  not  contain  sufficient  gradients  to 
force  adequate  adaption  to  the  natural  channel.  The  concentration  evident  in 
this  initial  grid  is  caused  by  the  boundary  point  distribution  through  the 
control  functions  in  the  original  WESCOR  formulation  as  mentioned  in  para¬ 
graph  6.  Another  example  is  the  grid  shown  in  Figure  9,  which  was  computed 
using  2-point  inverse  power  interpolation  with  the  exponent  set  to  be  10.  The 
grid  seems  to  be  trying  to  adapt  to  two  channels.  The  reason  is  not  apparent. 


Grid  Generation  Using  Hypothetical  Depths 

56.  The  next  effort  involved  using  an  adjusted  depth  field  in  which  the 
natural  channel  was  assumed  to  be  100  ft  deep,  with  zero  depths  specified  out 
of  the  channel.  In  addition,  as  Illustrated  in  Figure  10,  the  depths  were 
read  in  as  points  of  cross  sections  that  were  constructed  approximately  per¬ 
pendicular  to  the  center  line  of  the  channel.  A  total  of  approximately 

100  depth  points,  with  a  value  of  either  0  or  100,  were  input. 

57.  With  the  hypothetical  depth  field,  attraction  of  grid  points  to  the 
navigation  channel  was  achieved  for  most  of  the  interpolation  schemes.  The 
4-point  bilinear  with  smoothing  and  the  3-point  bilinear  with  no  smoothing 
schemes  resulted  in  unstable  computations .  However,  the  3-point  bilinear 
scheme  with  smoothing  was  stable  and  the  resulting  coordinate  system,  pre¬ 
sented  in  Figure  11,  contains  good  clustering  of  grid  points  in  the  channel. 
Figures  12  and  13  show  grids  generated  using  the  2-point  inverse  power  inter¬ 
polation  with  exponents  of  4  and  10,  respectively.  Little  impact  as  a  result 
of  the  change  in  the  exponent  is  observed.  The  grid  presented  in  Figure  14 
resulted  from  using  4-point  inverse  power  interpolation  with  an  exponent  of  4, 


while  that  in  Figure  15  used  a  Taylor  series  interpolation  with  2  points.  The 
latter  is  perhaps  the  most  physically  appealing  of  all  the  grids  generated. 

Influence  of  Weighting  Factors 

58.  All  of  the  grids  presented  up  to  now  were  generated  with  a  weight¬ 
ing  factor  of  1.0  for  each  of  the  three  contributing  integrals  in  Equation  10. 
To  demonstrate  the  relative  effect  of  each  factor,  Figures  16  and  17  are  pre¬ 
sented.  In  Figure  16,  the  concentration  weight  factor  has  been  set  to  0.0, 
with  the  smoothness  and  orthogonality  factors  retaining  a  value  of  1.0.  As 
expected,  no  attraction  to  the  navigation  channel  occurs  and  a  relatively 
smooth  grid  is  generated.  This  occurs  even  though  the  initial  grid,  shown  in 
Figure  8,  had  concentration  based  upon  boundary  point  distribution.  The 
variational  formulation  with  only  smoothness  and  orthogonality  smooth  out  the 
initial  concentration.  In  Figure  17.  both  the  concentration  factor  and  the 
orthogonality  factor  are  set  to  0.0,  while  keeping  the  smoothness  factor  at 
1.0.  This  is  the  grid  that  would  be  computed  by  the  WESC0R  code  with  zero 
values  for  the  control  functions  in  Equation  11.  An  interesting  observation 
is  that  the  grid  generated  with  no  enhancement  of  orthogonality,  i.e.,  Fig¬ 
ure  17,  would  probably  be  more  appropriate  for  flow  computations  than  the  one 
in  which  the  orthogonality  condition  is  enforced,  i.e.,  Figure  16,  because 
orthogonality  is  achieved  at  the  expense  of  smoothness. 

Practical  Aspects 

59.  Obviously,  not  only  the  numerical  grid,  in  which  grid  lines  follow 
navigation  channels,  but  also  the  water  depths  associated  with  each  computa¬ 
tional  cell  are  desired  as  output  from  the  grid  generator.  The  numerical 
hydrodynamic-transport  model  to  be  subsequently  employed  would  then  use  this 
output  in  the  computation  of  flow  and  constituent  fields.  However,  use  of  the 
actual  depth  field,  input  as  a  random  distribution,  did  not  result  in  suffi¬ 
cient  depth  gradients  to  force  grid  lines  to  adequately  follow  the  natural 
navigation  channel  in  Chesapeake  Bay.  As  a  result,  a  hypothetical  depth  field 
was  created,  with  large  depths  in  the  channel  and  zero  values  elsewhere.  In 
this  case,  the  final  depth  associated  with  each  computational  cell  is  meaning¬ 
less.  It  should  be  realized  that  as  far  as  grid  generation  is  concerned  the 
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use  of  either  real  or  adjusted  depths  is  immaterial,  although  real  depths  are 
required  on  the  final  grid  upon  which  flow  computations  are  to  be  made.  The 
solution  is  to  use  the  adjusted  depth  field  to  compute  the  grid  and  then  to 
interpolate  from  the  actual  depth  field  for  the  water  depth  to  be  associated 
with  each  cell  of  the  final  grid. 

60.  Based  upon  the  Chesapeake  Bay  grid  generation,  it  appears  that  the 
Taylor  series  interpolation  with  2  points  yields  perhaps  the  "best-looking" 
adaptive  grid.  Evidently  the  gradients  in  the  hypothetical  field  are  better 
maintained  after  interpolation  onto  the  initial  grid,  although  the  actual 
values  are  in  error.  Table  1  presents  the  depth  limits  resulting  from  dif¬ 
ferent  interpolation  schemes,  along  with  the  limits  of  the  interpolated  depth 
field  after  smoothing.  Based  upon  the  results  shown  in  Table  1,  it  would 
appear  that  interpolation  of  the  actual  depth  field  onto  the  final  adaptive 
grid  should  perhaps  be  based  upon  the  inverse  power  interpolat ion  with  no 
smoothing . 


PART  VI :  SUMMARY  AND  CONCLUSIONS 


Summary 

61.  As  a  result  of  the  need  for  generating  boundary-fitted  numerical 
grids  for  use  with  finite  difference  models  in  coastal  and  estuarine  areas,  a 
numerical  grid  generation  code  called  WESCORA  has  been  developed  as  an  exten¬ 
sion  of  an  earlier  code,  called  WESCOR.  Grid  generation  is  based  upon  adap¬ 
tive  grid  techniques  that  have  been  developed  for  dynamic  grid  adaption.  How¬ 
ever,  for  the  purpose  of  computing  free-surface  hydrodynamics  and  transport  in 
coastal  and  estuarine  areas  these  techniques  have  been  employed  here  only  for 
the  generation  of  the  initial  grid  and  no  time-dependence  is  considered. 

62.  Two  adaptive  mechanisms  are  contained  in  WESCORA;  namely,  adaption 
based  upon  a  variational  formulation  and  adaption  based  upon  a  control  func¬ 
tion  formulation.  In  both  cases,  adaption  to  water  depths  is  used  to  force 
grid  lines  to  follow  navigation  channels.  Preliminary  testing  of  the  two 
approaches  has  shown  that  the  variational  formulation  is  apparently  the  more 
reliable. 

63.  In  the  variational  approach,  a  functional  that  consists  of  the  sum 
of  three  integrals  involving  the  physical  coordinates  is  minimized  over  the 
grid.  The  first  of  the  integrals  controls  grid  point  concentration  while  the 
other  two  control  grid  smoothness  and  skewness  of  the  grid  lines.  By  weighing 
the  importance  of  the  three  integrals,  either  grid  point  concentration, 
smoothness,  or  orthogonality  can  be  emphasized. 

Conclusions 


64.  Preliminary  testing  of  WESCORm  on  Chesapeake  Bay  has  resulted  in 
the  general  conclusion  that  numerical  grids  can  be  generated  such  that  grid 
points  cluster  in  navigation  channels.  However,  unless  there  are  large  dif¬ 
ferences  between  the  channel  and  nonchannel  depths,  a  fictional  depth  field 
may  be  required  to  achieve  adequate  adaption.  In  this  case,  the  actual  depth 
field  should  be  interpolated  onto  the  final  adaptive  grid  to  provide  rela¬ 
tively  accurate  water  depths  to  be  associated  witn  each  computational  cell  cf 
the  grid.  Particular  conclusions  concerning  usage  of  the  code  are  offered 
below . 
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a_.  The  inverse  power  interpolation  is  bounded  by  the  extreme 
values  of  the  depth  values  and  is  the  smoothest  interpola¬ 
tion.  It  may  give,  however,  a  depth  distribution  that  is 
smoother  than  intended.  A  value  of  2  to  4  for  the  exponent 
involved  here  is  probably  appropriate. 

b_.  There  is  probably  no  real  reliable  advantage  to  using  large 

numbers  of  points  in  the  interpolation,  three  points  being  most 
reasonable . 

c_.  The  features  provided  requiring  a  convex  polygon  of  depth 

points  in  the  i nterpolation  and  smoothing  after  the  interpola¬ 
tion  should  be  used. 

d_.  Multiple  groups  of  interpolation  points  should  be  used  with  the 
Taylor  series  and  bilinear  interpolations,  four  groups  being 
reasonable . 

_e .  For  grid  generation,  the  bilinear  interpolation  with  three 
points  or  the  Tavlor  series  interpolation  with  two  points, 
should  be  good  choices  in  most  cases.  The  inverse  power'  in¬ 
terpolation  with  no  smoothing  should  provide  relatively 
accurate  depths  on  the  final  grid. 

f_.  tJeumann  boundary  conditions  (orthogonality  at  the  boundary) 
should  be  used  on  boundary  segments  on  which  the  depth 
varies.  This  feature  requires  that  the  boundary  segment  in 
question  be  straight. 

g_.  The  variational  .adaptive  mechanism  is  the  more  reliable.  The 
values  of  the  smoothing  and  orthogonality  parameters  should  be 
kept  equal  and  not  more  than  an  order  of  magnitude  smaller  than 
the  concentration  parameter  in  most  cases. 

h_.  If  the  randomly  distributed  depth  field  does  not  yield  depth 
gradients  sufficient  to  cause  good  adaption  to  navigation 
channels,  a  hypothetical  field  such  as  the  one  used  in  the 
Chesapeake  Bay  example  should  be  input.  In  general,  cross 
sections  should  be  aligned  perpendicular  to  the  channel. 
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Figure  3-  Inverse  power  interpolation 
among  depth  points 
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Figure  13.  Depth^adapti ve  grid  u3ing  inverse  power  interpolation 
with  two  points  and  an  exponent  of  10 


Figure  14.  Depth-'-adapti  ve  grid  using  inverse  power  interpolation 
with  four  points  and  an  exponent  of  2 


Figure  15.  Depth-*adapti ve  grid  using  Taylor  series  interpolation 

with  two  points 


Figure  16.  Depth-adaptive  grid  with  concentration  factor  =  OJ 
and  orthogonality  and  smoothness  factors  =  1.0 


Figure  17.  Depth-adaptive  grid  with  concentration  and 
orthogonality  factors  =  0.0  and  smoothness  factor  =  1.0 


Table  1 

Interpolated  Depth  Limits 


Interpolation  Scheme 

Actual  Limits 

Bilinear  with  3  points 

0-156 

Bilinear  with  3  points 

0-100 

Inverse  power  with 

4  points,  exponent  = 

4 

0-156 

Inverse  power  with 

4  points,  exponent  = 

2 

0-100 

Inverse  power  with 

3  points,  exponent  = 

4 

0-156 

Inverse  power  with 

3  points,  exponent  = 

4 

0-100 

Taylor  series  with 

2  points 

0-156 

Taylor  series  with 

2  points 

0-100 

Taylor  series  with 

3  points 

0-156 

Interpolated  Limits  Smoothed  Limits 


0-120 

0-  86 

0-100 

0-  98 

0-137 

0-100 

0-100 

0-  95 

0-137 

0-100 

0-100 

0-100 

0-139 

1-  96 

0-449 

1-296 

0-204 

1-126 

APPENDIX  A:  INPUT  INSTRUCTIONS 
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C 

cttmtniiitmmmttK  u  e  s  c  o 


2-0  BOUNDARY-FITTED  COORDINATE  SYSTEM  CODE  (  DEPTH-ADAPTIVE  > 
MISSISSIPPI  STATE  UNIOERSITV  .  19«2  REVISED  1985 


U.S.  ARMY  ENGINEER  WATERWAYS  EXPERIMENT  STATION 
VICKSBURG,  MISSISSIPPI 


cmtmimitttiititiiiiiitititiimtimtimtmmtmtmmttimt 

c 

C  tttmmttll  SLIT-SIAI  CONFIGURATION  ttttt 

s 

tm  ATTRACTION  TO  COORDINATE  LINES/POINTS  AND  TO  SPACE  LINES/POINTS. 

»tt  CONTROL  FUNCTIONS  ALSO  INTERPOLATED  FROM  IOUNDARY  POINT  DISTRIBUTION 


tttt  ADAPTIVE  ON  DEPTH  DISTRIBUTION, 
t 

tttt  CAN  BE  ORTHOGONAL  ON  STRAIGHT  BOUNDARY  SEGMENTS. 


ttttttxsttttxtttttxttttttxtttttttttxttttttttttxxtttxtttttttttttttttt 
ttttttttttt  INPUT  INSTRUCTIONS  * 
t 

ttt  CARDS ( 2 )  <  LABEL  -  FORMAT (18A8) 
t 

t  LABEL  -  TWO  8B  CHARACTER  CARDS.  (BLANK  CARDS  IF  NO  LABEL) 


ttt  CARD 
t 


IMAX. JMAX.NBDV. ITER. ISLIT. IBNDRV, IDISK, IUIR, IWINTL, 
IUFIN.NREN.NEUINT  -  FORMATC 1215 ) 


t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 


IMAX  -  NUMBER  OF  XI  POINTS. 


JMAX  -  NUMBER  OF  ETA  POINTS. 


NBDY  -  TOTAL  NUMBER  OF  SLAB  SIDES  AND  SLITS  IN  THE  FIELD. 


ITER  -  MAXIMUM  NUMBER  OF  ITERATIONS  ALLOWED. 


ISLIT  -  -1  SLAB  SIDES  OR  SLITS  READ  FROM  CARDS. 

X.V  -  FORMAT <2F IB. B)  ,  ONE  POINT  PER  CARD. 
*2  SLAB  SIDES  OR  SLITS  READ  FROM  FILE  19. 

X.V  -  FORMAT( 2F19.B )  ,  ONE  POINT  PER  CARD. 


(NOTE*  HORIZONTAL  SLITS  ARE  READ  CLOCKWISE  FROM  RIGHT  END.) 


VERTICAL  SLITS  ARE  COUNTER-CLOCKWISE  FROM  TCP. 
SLAB  SIDES  MAY  BE  READ  IN  EITHER  DIRECTION. 


IBNDRV  - 


•0  OUTER  BOUNDARY  CALCULATED  INTERNALLY  AS  CIRCLE. 
•1  OUTER  BOUNDARY  READ  FROM  CARDS. 

-  FORMAT < 2F 19 .9 )  ,  ONE  POINT  PER  CARD. 


X.  Y 


•2  OUTER  BOUNDARY  READ  FROM  FILE  19. 


X.Y 


F0RMAT(2F19.9)  .  ONE  POINT  PER  CARD. 


— 1  OUTER  BOUNDARY  READ  IN  SEGHENTS  AS  SLAB  SIDES. 


(NOTE*  FOR  IBNDRV  •  1  OR  2  ,  OU^ER  BOUNDARY  IS  READ  CLCCKUISE 


FROM  POINT  ( INFXI , INFETA ) . 


i 


i 


i 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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•OUTER  BOUNDARY •  REARS  ENTIRE  BOUNDARY  OF  TRANSFORMED 
REGION  IF  NREN-0.  IF  NREN  IS  NO*  ZERO,  THEN  OUTER 
BOUNDARY  IS  THE  TOR  OF  THE  TRANSFORMED  REGION  AND 
INNER  BOUNDARY  IS  THE  BOTTOM. 


IDISK 


-  ••  DON'T  READ  OR  WRITE  SYSTEM  FROM  OR  ON  FILE. 

•1  WRITE  SYSTEM  Oh  FILES  10  I  11.  DON'T  READ  SYSTEM  FR 
•2  WRITE  SYSTEM  ON  FILES  IB  1  11  .  READ  SYSTEM  FROM  FI 
•3  READ  SYSTEM  FROM  FILE  10  FOR  RESTART.  DON'T  WRITE  S 


(NOTES  FILE  10  IS  RESTART  FILE  FOR  CONTNUATION  OF  ITERATION.) 
(  FILE  11  IS  STORACE  FILE  FOR  FINAL  SYSTEM.  ) 


IWIR 


-  -0  DON'T  PRINT  EACH  ITERATION  ERROR. 
•I  PRINT  EACH  ITERATION  ERROR. 


IWINTL  -  *0  DON'T  PRINT  INITIAL  GUESS. 
•1  PRINT  INITIAL  GUESS. 


IUFIN 

NREN 


-  ZERO  SUPPRESSES  PRINT  OF  FINAL  VALUES. 


NON-ZERO  USES  RE-ENTRANT  BOUNDARY  ON  LEFT  l  RIGHT  SIDE 
OF  TRANSFORMED  REGION,  WITH  OUTER  BOUNDARY  ON  TOP 
AND  INNER  BOUNDARY  ON  BOTTOM. 


INNER  BOUNDARY  IS  READ  AS  FOLLOWS  BEFORE  READING  OUTER 
•1  INNER  BOUNDARY  READ  FROM  CARDS. 

X,V  -  FORMAT(2F10.0)  ,  ONE  POINT  PER  CARD. 

-2  INNER  BOUNDARY  READ  FROM  FILE  10. 

X,Y  -  FORMAT (2F 10.0)  ,  ONE  POINT  PER  CARD. 


(NOTES  SLITS  AND/OR  SLABS  MAY  ALSO  BE  PRESENT.) 


NEUINT  -  1  SUPPRESSES  NEUMANN  BOUNDARY  CONDITIONS  FOR 
INITIAL  GRID. 


M  CARDS(NBDY)  LB1,LB2.LB3.LTYPE,NCU»0D 


FORMAT(SIS) 


LB1.LB2  -  FIRST  AND  LAST  INDICES  OF  SLAB  SIDE  OR  SLIT  ENDS. 

(LB2  NAY  BE  LESS  THAN  LB1  FOR  SLAB  SIDE.  INPUT  IS  FRO 


LB3  -  INDEX  OF  LINE  ON  WHICH  SLAB  SIDE  OR  SLIT  IS  LOCATED. 
LTVPE  -  SLAB  SIDE  OR  SLIT  TYPE  (1  FOR  HORIZONTAL,  2  FOR  UERTICA 


(NEGATIVE  INDICATES  SLAB  SIDE,  RATHER  THAN  SLIT.) 
(SUBTRACT  10  FOR  OUTER  BOUNDARY  SEGMENT.  ) 
(I.E.,  -11  IS  HORIZONTAL  OUTER  BOUNDARY  SEGMENT,) 
(  -12  IS  VERTICAL  OUTER  BOUNDARY  SEGMENT.  ) 


NEUBOD  -  1  MEANS  NEUMANN  BOUNDARY  CONDITION. 
0  MEANS  DIRICHLET  CONDITION. 


SS  CARO 


NEUOUTM)  -  FORMAT (4 15  ) 


NEUOUT  -  1  MEANS  NEUMANN  BOUNDARY  CONDITIONS  ON  OUTER  BOUNDARY. 
(OUTER  BOUNDARY  MADE  OF  SLAB  SIDES  IS  CONTROLLED 
BY  NEUBOD  INSTEAD  OF  NEUOUT.) 

0  MEANS  DIRICHLET  CONDITIONS. 
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(1)  »  LEFT 

(2)  «  RIGHT 

(3)  t  BOTTOM 

(4)  t  TOP 

St  CARD  >  NDEP,NDOUT,NDOP,NDEPF.NlNP,NCOM,NTES,NSMO  -  F0RMAT(8IS> 

NDEP  -  NUMBER  OF  DEPTH  POINTS. 

(  ZERO  FOR  NON- ADAPTIVE  SYSTEM  ) 

NDOUT  -  NON-ZERO  ADDS  OUTER  BOUNDARv  POINTS  TO  DEPTH  POINT 
LIST,  UITH  ZERO  DEPTH. 

<  THIS  SHOULD  BE  USED  UNLESS  THE  INPUT  DEPTH  POINTS 
EXTEND  TO  OR  BEVOND  THE  OUTER  BOUNDARY.  INPUT  DEPTH 
POINTS  THAT  ARE  COINCIDENT  WITH  OUTER  BOUNDARV  POINTS 
WILL  RETAIN  THE  INPUT  DEPTH  EUEN  WHEN  THIS  FEATURE 
IS  USED.  ) 

NDOP  -  1  PRINTS  DEPTH  POINTS. 

NOEPF  -  NON-ZERO  READS  DEPTH  POINTS  FROM  FILE  12. 

NINP  -  INTERPOLATION  TYPE  *  9  >  HARMONIC  -  ALL 

13  >  BILINEAR  -  3  POINTS 

14  >  BILINEAR  -  4  POINTS 

2  >  2-POINT  TAYLOR  SERIES 

3  >  3-POINT  TAYLOR  SERIES 

4  >  4-POINT  TAYLOR  SERIES 

5  >  S-POINT  TAYLOR  SERIES 
-2  )  HARMONIC  -  2  POINTS 
-3  >  HARMONIC  -  3  POINTS 
-4  >  HARMONIC  -  4  POINTS 
-S  >  HARMONIC  -  5  POINTS 

l  BEST  ARE  13  ,  2  ,  3  ,  -2  ,  -3  ) 

NCON  -  NUMBER  OF  INTERPOLATION  POINT  SETS  FOR  TAYLOR.  <  4  ) 
POWER  FOR  HARMONIC.  (  2  ) 

NTES  -  NON-ZERO  REQUIRES  SUCCEEDING  POINTS  TO  SURROUND 
INTERPOLATION  POINT.  (  USE  THIS  ) 

HSMO  -  1  SUPPRESSES  DEPTH  SMOOTHING.  C  DON'T  SUPPRESS  ) 

St  CARD  R<1 >.R(2),R(3),YINFIN,AINFIH,X0INF.Y0INF,INFXI.INFETA 
-  FORMAT ( 7F 1 6 . f , 2 I 5 ) 

R( 1 )  -  SOR  ACCELERATION  PARAMETER  FOR  INITIAL  GRID. 

(ZERO  VALUE  CAUSES  UARIABIE  ACCELERATION  PARAMETER) 
(FIELD  TO  BE  CALCULATED  INFERNALLY .  > 

(NEGATIUE  VALUE  GIVES  ALGEBRAIC  SYSTEM) 

R<  2  )  -  ALLOWABLE  X  ITERATION  ERROR. 

R( 3  )  -  ALLOWABLE  V  ITERATION  ERROR. 

YINFIN  -  RADIUS  OF  CIRCULAR  OUTER  BOUNDARV. 

AINFIN  -  ANGLE  OF  FIRST  POINT  ON  CIRCULAR  OUTER  BOUNDARV  (DEGRE 


ooooooooooooooooooooooonoooooooooooooooonoooooooooooooooooo 
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(COUNTER-CLOCK  FROM  POSITIVE  X-AXIS.) 

XOINF.YOINF  -  CENTER  OF  CIRCULAR  OUTER  10UNDARV. 

INFXI.INFETA  -  INDICES  OF  FIRST  POINT  ON  CIRCULAR  OUTER  BOUNDAR 

(NOTE  t  LAST  7  OF  THESE  PARARETERS  ARE  IRRELEVANT  IF  OUTER  B 

X X  CARD  »  ACC L.WF AC  I, SFAC, OF AC  I  -  F ORMAT ( 4F1B . 0  ) 

(  OMIT  IF  NDEP»B  ) 

ACCL  -  SOR  ACCELERATION  PARANETER  FOR  ADAPTIVE  CRID.  (  1.8  > 

UF AC  I  -  CONCENTRATION  FACTOR.  (  B.0-1.0  > 

(  POSITIVE  VALUE  FOR  VARIATIONAL  ADAPTION  ) 

(  NEGATIVE  VALUE  FOR  CONTROL  FUNCTION  ADAPTION  ) 

(  SFAC  1  OFACI  ARE  IRRELEVANT  IN  UFACI  IS  NEGATIVE) 

SFAC  -  SMOOTHNESS  FACTOR.  (  B.B-1.8  ) 

OFACI  -  ORTHOGONAL  I  TV  FACTOR.  (  B.B-1.0  ) 


X  IF  BODIES  AND/OR  OUTER  BOUNDARY  ARE  READ  FRON  CARDS.  SUCH  CARDS 
X  FOLLOW  NEXT  UNLESS  A  RESTART  IS  USED. 

S 

X  SLITS  AND/OR  SLAB  SIDES  ARE  READ  FIRST,  THIN  OUTER  BOUNDARY  IS  READ. 
(THESE  RULES  APPLY  FOR  READING  FRON  FILE  IB  AS  WELL  AS  FRON  CARDS.) 

X 


X  IF  DEPTH  POINTS  ARE  READ  FRON  CARDS,  SUCH  CARDS  FOLLOU  NEXT 
X  IN  THE  FOLLOWING  FORMAT.  IF  DEP-B  NO  DEPTH  POINTS  ARE  READ. 

X  ONIT  IF  A  RESTART  IS  USED,  UNLESS  IT  IS  FRON  A  CONVERGED  INITIAL 
X  GRID  GENERATED  WITHOUT  READING  DEPTH  POINTS. 


XX  CARDS (NDEP  )  I  XDEP. YDEP, DEP  -  FORMAT ( 3F1B.8 > 

XDEP.VDEP  -  X , Y  COORDINATES  CF  DEPTH  POINT. 
DEP  -  DEPTH. 


:X  IF  NO  COORDINATE  ATTRACTION  IS  TO  BE  USED,  FOLLOW  THESE  CARDS 
X  WITH  FIVE  BLANK  CARDS.  IF  ATTRACTION  IS  TO  BE  USED.  USE  THE  FOLLOWING 
X  INPUT  RATHER  THAN  THE  BLANK  CARDS: 

X 

X  INPUT  FOR  COORDINATE  SYSTEM  CONTROL  :  USE  FOUR  SETS,  ONE  fOR 
X  XI -LINE  ATTRACTION  ^0  COORDINATE  LINES/POINTS,  ONE  FOR  ETA-LINE  ATTRAC 
X  TO  COORDINATE  LINES/POINTS,  ONE  FOR  XI -LINE  ATTRACTION  to  SPACE  LINES/ 
X  AND  ONE  FOR  ETA-LINE  ATTRACTION  TO  SPACE  LINES/POINTS. 

X  ANY  SET  NOT  WANTED  IS  REPLACED  BY  ONE  BLANK  CARD. 

X 


c  SB 

C  *t  THE  FOLLOWING,  NARKED  WITH  #.  IS  FOR  ATTRACTION  TO  COORDINATE  LlN£S/PO 

C  It 

C  StSS  CARO  »  ATVP, ITVP, NLN, NPT, DEC, ANPFAC  -  FORNAT( A8, IS, 215, SF 19. A) 

C  St 

C  St  AT VP  -  TYPE  OF  ATTRACTION.  (XI  FOR  XI-LINE  ATTRACTION, 

C  St  ETA  FOR  ETA-LINE  ATTRACTION.)  LEFT  JUSTIFIED. 

C  St 

C  St  ITVP  -  ZERO  GIUES  ATTRACTION  ON  BOTH  SIDES. 

C  St  NON-ZERO  GIUES  ATTRACTION  ON  UPPER  SIDE  AND 

C  St  REPULSION  ON  LOWER  SIDE. 

C  St 

C  St  NLN  -  NUNBER  OF  ATTRACTION  LINES. 

C  St 

c  St  NPT  -  NUNBER  OF  ATTRACTION  POINTS. 

C  St 

c  St  DEC  -  NON-ZERO  DEC  USES  DEC  FOR  DECAV  FACTOR. 

c  St 

C  St  ANPFAC  -  NON-ZERO  ANPFAC  NULTIPLIES  ALL  AMPLITUDES  BY  ANPFAC. 

C  St 

C  StSS  CARDS (NLN )  :  JLN, ALN, DLN  -  F0RNAT( SX, IS, 2Flt. 0 ) 

C  St  (ONIT  IF  NLN  IS  ZERO) 

C  St 

C  St  JLN  -  ATTRACTION  LINE  INDEX. 

C  St 

C  St  ALN  -  ANPLITUDE  (NEGATIVE  REPELS)  FOR  LINE  ATTRACTION. 

C  St 

C  St  DLN  -  DECAV  FACTOR  FOR  LINE  ATTRACTION. 

C  St 

C  StSS  CARDS ( NPT )  »  IPT, JPT.APT,DPT  -  FORNAT(2IS,2F10.8 ) 

C  St  (ONIT  IF  NPT  IS  ZERO) 

C  St 

C  St  IPT.JPT  -  ATTRACTION  POINT  INDICES. 

C  St 

C  St  APT  -  ANPLITUDE  (NEGATIVE  REPELS)  FOR  POINT  ATTRACTION. 

C  St 

C  St  DPT  -  DECAV  FACTOR  FOR  POINT  ATTRACTION. 

C  ss 

C  Sftffttftftfttittttttttfttttttlttftfttftttttftttttttttttttttttttttt 
C  St 

C  Stt  THE  FOLLOWING,  NARKED  WITH  f,  IS  FOR  ATTRACTION  TO  SPACE  LINES/POINTS 

C  St 

C  StSS  THE  FOLLOWING  CARDS  ARE  FOR  ATTRACTION  TO  LINES  AND/OR  POIN-S 
C  StSS  DEFINED  BY  X,V  COORDINATES.  IF  NLN  IS  NOT  ZERO,  THEN  NLN 

C  StSS  OF  THE  CARDS  GIVING  NP  NUST  APPEAR.  EACH  OF  THESE  CARDS  IS 

C  StSS  FOLLOWED  BY  NP  OF  THE  CARDS  GIVING  XPT,  ETC.  IF  NPT  IS  NOT 

C  StSS  ZERO,  THEN  NPT  OF  THE  CARDS  GIVING  XPT.  ETC.  NUST  FOLLOW 

C  StSS  THE  LAST  GROUP  OF  THESE  CARDS. 

C  StSS  ANY  SET  NOT  WANTED  IS  REPLACED  BY  ONE  BLANK  CARD. 

C  St 

C  StSS  CARD  l  ATYP, ITVP, NLN, NPT, DEC, ANPFAC  -  FORNAT ( A8. 12. SIS. 2F10 .9 > 


c 

St 

c 

St 

ATVP 

-  TYPE  OF  ATTRACTION.  (XI  FOR  XI-LINE  ATTRACTION, 

c 

St 

ETA  FOR  ETA-LINE  ATTRACTION.)  LEFT  JUSTIFIED. 

c 

St 

c 

St 

ITYP 

-  ZERO  GIUES  ATTRACTION  OH  BOTH  SIDES. 

c 

St 

NON-ZERO  GIVES  ATTRACTION  ON  UPPER  SIDE  AND 

c 

St 

REPULSION  ON  LOWER  SIDE. 

c 

St 

A  7 
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NLN  -  NUMBER  Of  ATTRACTION  tlN£S. 

NPT  -  NUMBER  Of  ATTRACTION  POINTS. 

(NOT  INCLUDING  POINTS  ON  ATTRACTION  LINES) 

DEC  -  NON-ZERO  DEC  USES  DEC  fOR  DECAY  FACTOR. 

AMPFAC  -  NON-ZERO  AMPFAC  MULTIPLIES  ALL  AMPLITUDES  1Y  AMPFAC. 


tt  CARD  >  HP 


FORMAT ( 15  ) 


NP  -  NUMBER  Of  POINTS  ON  THIS  ATTRACTION  LINE. 


tt  CARDS  :  XPT,YPT,APT.DPT,UEC1,U£C2 


FORMAT ( 6P10 . 0  ) 


XPT , YPT  -  COORDINATES  OF  ATTRACTION  POINT  OR 
POINT  ON  ATTRACTION  LINE. 

APT  -  ATTRACTION  AMPLITUDE  (NEGATIUE  REPELS). 

DPT  -  DECAY  FACTOR. 

VEC1.VEC2  -  X , Y  COMPONENTS  OF  UNIT  UECTOR  NORMAL  TO 

ATTRACTION  DIRECTION  FOR  POINT  ATTRACTION. 
(CALCULATED  INTERNALLY  FOR  LINE  ATTRACTION.) 


tt  THE  LAST  COORDINATE  SYSTEM  CONTROL  CARD  IS  THE  FOLLOWING  CARD  « 

It  CARD  :  IFAC. IRIT.EFAC  -  FORMAT (2IS, F 10. 0  ) 

( CAN  BE  USED  TO  AID  CONVERGENCE  BV  CONVERGING  FIELD  ) 

(WITH  LESS  ATTRACTION  FIRST  AND  USING  THIS  RESULT  ) 

(AS  THE  INITIAL  GUESS  FOR  STRONGER  ATTRACTION.  ) 

(BLANK  CARD  MUST  BE  INPUT  IF  THIS  FEATURE  IS  NOT  USED.) 
(STANDARD  IS  TO  NOT  USE  THIS  FEATURE  ,  BUT  ITS  USE  MAY) 

(BE  NECESSARY  UITH  STRONG  ATTRACTION.  ) 

IFAC  -  NUMBER  Of  STEPS  IN  ADDITION  OF  INHOMOGENEOUS  TERM. 
DOUBLES  INHOMOGENEOUS  TERM  AT  EACH  STEP. 

(ZERO  CONVERGES  WITH  FULL  ATTRACTION. 

(1.0  CONVERGES  UITH  NO  ATTRACTION  FIRST,  THEN  ) 

(UITH  FULL  ATTRACTION.  2.0  CONVERGES  WITH  NO  ) 

(ATTRACTION  FIRST,  THEN  UITH  HALF,  THEN  UITH  FULL.) 
(INCREASE  NUMBER  OF  STEPS  IF  DIVERGENCE  OCCURS.  ) 

IRIT  -  NON-ZERO  VALUE  CAUSES  INHOMOGENEOUS  TERM  TO  BE  PRINTED. 

EFAC  -  MULTIPLE  OF  CONVERGENCE  CRITERION  *0  BE  USED  FOR 
INTERMEDIATE  CONVERGENCE  BETWEEN  ADDITIONS  OF 
INHOMOGENEOUS  TERM.  (TYPICALLY  10.0  ) 

It  CARD  »  LDEP,LDOUT,LDOP,LDEPF,LINP,LCCM,LTES.LSMO  -  FORMAT ( 815  > 

(  THIS  CARD  ALLOUS  A  SEPARATE  SET  OP  DEPTH  POINTS  ~0  BE  INPUT 
TO  BE  USED  FOR  INTERPOLATION  OF  DEPTHS  ON  THE  PINAL  GRID  '0 


*  •  "  »  *  «  *  .  *  «  *  ■  •  ”  •  *  -  T.  ***»“*  •-*_  .  *  ,  •  .  .  *  *  ' 
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*  t€  OUTPUT  ALONG  UITH  THE  GRID.  ^1$  IS  USEFUL  in  THE  CASE 

X  WHERE  ADJUSTED  DEPTHS  WERE  USED  FOR  THE  ADAPTION  WHILE  ACTUAL 

t  DEPTHS  ARE  DESIRED  FOR  USE  UITH  THE  GRID  IN  F LOU  CODES. 

*  INPUT  A  BLANK  CARD  HERE  IF  THE  FIRST  SET  OF  DEPTH  POINTS  IS 

X  TO  BE  USED  FOR  THE  OUTPUT  DEPTHS.  IN  EITHER  CASE  THE  OUTPUT 

*  DEPTHS  ARE  LOCATED  AT  CELL  CENTERS  BV  AUERAGING  THE  FOUR 

X  SURROUNDING  GRID  POINT  UALUES.  THE  <I.J>  OUTPUT  DEPTH  IS  FOR 

X  THE  CELL  WITH  <I,J)  AT  THE  LOWER-LEFT  CORNER.  ) 

*  LDEP  -  NUNBER  OF  DEPTH  POINTS. 

X  t  ZERO  FOR  NON-ADAPTIUE  SYSTEN  > 

X  LDOUT  -  NON-ZERO  ADDS  OUTER  BOUNDARY  POINTS  TO  DEPTH  POINT 

X  LIST,  WITH  ZERO  DEPTH. 

t  (  THIS  SHOULD  BE  USED  UNLESS  THE  INPUT  DEPTH  POINTS 

t  EXTEND  TO  OR  BEYOND  THE  OUTER  BOUNDARY.  INPUT  DEPTH 

X  POINTS  THAT  ARE  COINCIDENT  WITH  OUTER  BOUNDARY  POINTS 

X  WILL  RETAIN  THE  INPUT  DEPTH  EUEN  WHEN  THIS  FEATURE 

X  IS  USED.  ) 

X  LOOP  -  1  PRINTS  DEPTH  POINTS. 

X  LDEPF  -  NON-ZERO  READS  DEPTH  POINTS  FROfl  FILE  16. 

X  LINP  -  INTERPOLATION  TVPE  i  9  >  HARNONIC  -  ALL 

X  13  >  BILINEAR  -  3  POINTS 

X  14  >  BILINEAR  -  4  POINTS 

X  2  >  2-POINT  TAYLOR  SERIES 

X  3  >  3-POINT  TAYLOR  SERIES 

X  4  >  4-POINT  TAYLOR  SERIES 

X  5  >  S-POINT  TAYLCR  SERIES 

X  -2  >  HARNONIC  -  2  POINTS 

X  -3  >  HARNONIC  -  3  POIN’S 

X  -4  >  HARNONIC  -  4  POINTS 

X  -5  >  HARNONIC  -  S  POINTS 

X  <  BEST  ARE  13  ,  2  .  3  ,  -2  ,  -3  ) 

X  LCON  -  NUNBER  OF  INTERPOLATION  POINT  SETS  FOR  TAYLOR.  <  4  ) 

X  POWER  FOR  HARNONIC.  (  2  ) 

X  LTES  -  NON-ZERO  REQUIRES  SUCCEEDING  POINTS  TO  SURROUND 

X  INTERPOLATION  POINT.  (  USE  THIS  ) 

X  LSHO  -  l  SUPPRESSES  DEPTH  SNOOTHING. 

XXX  CARDS ( LDEP  >  :  XDEP, VDEP, DEP  -  FORNATl 3F19.0  ) 

X  XDEP.YDEP  -  X.Y  COORDINATES  OF  DEP’X  POINT. 

X  DEP  -  DEPTH. 

XtXXXXXXXXXIXXXXXlXXXXXXXXXXXXXXXXXlXtXXSXXlXXXXtStSXSXSttXSXXXXXfXXX 
X  NASS  STORAGE  FILES  > 

*  RESTART  FILE  -  READ  FRON  FILE  14  ,  WRITTEN  ON  FILE  15  • 

*  BFSTART  CAN  BE  FRON  PARTIALLY  OR  FULLY  CONUERCED 

X  INITIAL  OR  ADAPTIUE  GRID. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


(  14)  RXl.RETA.J.uL.WU.RFAC.OfAC.DEP.XDEP.VDEP 
(14)  X.V.LSLIT. LABEL, I  RAX, JRAX 

(  14  )  NBDV. HUM.  LB1»LB2.LB3.LTVP£,  LPT.  XL,  XU,  VL.  YU, 

DIAI.DIRJ.DIIW.DIRB. LACC.UACC.NEUOUT.NEUBOD, I T 1 

COORDINATE  SVSTEA  STORAGE  FILE  -  FILE  11  i 

(11)  LABEL,  I1AX,  JflAX 

(11)  (  (  LSLIT  (  I ,  J  ) ,  I  •  1  ,  INAX  ) ,  J  *  1 .  JflAX  ) 

(11  )  ((X(I,J),I«1,I  F1AX  ) ,  J  »  1 ,  JflAX  ) 

(11  )  (  (V(I,J),I«1, IflAX  ) , J« 1 , JflAX  ) 

( 11  )  (  <U(  I,  J  >,1-1, 1  WAX  ),  J-l  .  JflAX  ) 

(OflITTED  IF  NDEP*B) 

(11 )  NBDV , MUNI, LB1 , LB2 , LB3, LTVPE , LPT , XL, XU, VL, VU, 
DIRI,DIHJ,DIflP.DIAB 
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The  following  files  are  defined: 


BAYR1 

BAYD1 

BAY1 

RPLOTA 

CSPLOT 

WESCORA 


CDC  run  stream  file. 

File  containing  (x,y)  location  and  the  value  of  depth  points. 
One  depth  point  per  line. 

File  containing  (x,y)  location  of  boundary  points.  One  point 
per  line. 

CDC  plot  run  stream. 

Plot  code  source. 

Grid  code  source. 


BAYR1 

/JOB 

JFTlXXX(Tf  .CH3MB—.M) 
USER,  .  • KOC ■ 

CHARGE, 

SBULIH.AB-7????. 

GET, OWCSCOR. 
GET.TAPE18-BAVDI. 

GET. TAPE t6*»AV08. 

GET, TAPE1#*MV1 . 

OWCSCOR. 

REPLACE , OUTPUT • OUT . 

REPLACE . TAPE 1 5 • RE START . 
REPLACE. TAPE11*S0LN. 

DAVE ILE. DAY. 

REPLACE, DAY. 

EXIT. 

REPLACE, OUTPUT -OUT. 

DAYF I LE ,  DAY . 

REPLACE. DAY. 

/EOR 

ADAPTIVE 
CHESAPEAKE  BAY 


BAY1 


/s. 


RPLOTA 


'JOB 

8HJtXXX(CH37770e,Tl#,P4 ) 

USER,  ,  ,KOE . 

charge 

GET ,GCS, 0LD*0CSPL0T , OPLTFIL . 
COPYL, OLD, OPLTFIL, OCSPLOT, . A9 . 
REyiND.OCSPLOT . 

GET , TAPE  1 0»SOLN . 

RFL,  14IM0. 

CALL.GCSU-0CSPL0T,DEy.QCST4l  > 
REPLACE. TAPES6*PL0TCS. 
DAVFILE.DAY. 

REPLACE. DAY. 

EXIT. 

REPLACE . OUTPUT -BADOUT . 
DAVFILE.DAY. 

REPLACE. DAV. 


As  can  be  seen,  the  day  file  from  both  a  grid  generation  run  as  well  as 
a  run  to  generate  a  coordinate  plot  is  stored  in  a  permanent  file  called  DAY, 
Printed  output  from  a  grid  run  is  stored  in  a  permanent  file  called  OUT.  The 
following  sequence  of  commands  would  be  issued  from  a  Tektronix  terminal  to 
first  run  the  grid  model,  then  run  the  plot  code,  and  finally  to  plot  the 
boundary-fitted  coordinate  system  on  the  screen. 

Creates  a  local  file  called  BAYR1  from  the  permanent 
file  BAYR1 . 

Run  stream  BAYR1  is  submitted  as  a  batch  job. 

Lock  at  day  file  to  see  if  the  grid  run  was  successful. 

Look  at  the  printed  output  from  the  grid  run. 

Creates  a  local  file  called  RPLOTA  from  the  permanent 
file  RPLOTA. 

Run  stream  RPLOTA  is  submitted  as  a  batch  job. 

Look  at  day  file  to  see  if  the  plot  run  was  successful. 

Coordinate  system  is  plotted  on  the  screen 


GET, BAYR1  - 

SUBMIT, BAYR1  - 
XEDIT , DAY , P  - 
XEDIT , OUT , P  - 
GET, RPLOTA  - 

SUBMIT, RPLOTA  - 

XEDIT, DAY, P  - 

OLD, PLOTCS  | 
LHN  I 


